arXiv:1508.07583vl [hep-lat] 30 Aug 2015 


INT-PUB-15-042, NSF-KITP-15-120, NT@UW-15-09, MIT-CTP-4709 


Two Nucleon Systems at ^ 450 MeV from Lattice QCD 

Kostas Orginos/’^ Assumpta Parreno,^ Martin J. Savage,^ 

Silas R. Beane,^ Emmanuel Chang,^ and William Detmold^ 

(NPLQCD Collaboration) 

^Department of Physies, College of William and Mary, Williamsburg, VA 23187-8795, USA 
^Jefferson Laboratory, 12000 Jefferson Avenue, Newport News, VA 23606, USA 
^Dept. dDstruetura i Constituents de la Materia. Institut de Cieneies del Cosmos (ICC), 
Universitat de Bareelona, Marti Franques 1, E08028-Spain 
Institute for Nuelear Theory, University of Washington, Seattle, WA 98195-1550, USA 
^Department of Physies, University of Washington, Box 351560, Seattle, WA 98195, USA 
^ Center for Theoretieal Physies, Massaehusetts Institute of Teehnology, Cambridge, MA 02139, USA 

(Dated: September 1, 2015) 

Nucleon-nucleon systems are studied with lattice quantum chromodynamics at a pion 
mass of ^ 450 MeV in three spatial volumes using n/ = 2 + 1 flavors of light quarks. 

At the quark masses employed in this work, the deuteron binding energy is calculated to 
be Bd = 14.412 6 MeV, while the dineutron is bound by B^n = I2.5I5 0 MeV. Over the 
range of energies that are studied, the S-wave scattering phase shifts calculated in the ^Sq 
and ^Si-^Di channels are found to be similar to those in nature, and indicate repulsive short- 
range components of the interactions, consistent with phenomenological nucleon-nucleon 
interactions. In both channels, the phase shifts are determined at three energies that lie 
within the radius of convergence of the effective range expansion, allowing for constraints to 
be placed on the inverse scattering lengths and effective ranges. The extracted phase shifts 
allow for matching to nuclear effective field theories, from which low energy counterterms 
are extracted and issues of convergence are investigated. As part of the analysis, a detailed 
investigation of the single hadron sector is performed, enabling a precise determination of 
the violation of the Gell-Mann-Okubo mass relation. 

PACS numbers: 11.15.Ha, 12.38.Gc, 


I. INTRODUCTION 

Calculating the interactions between nucleons and the properties of multi-nucleon systems di¬ 
rectly from quantum chromodynamics (QCD) will be an important milestone in the development 
of nuclear physics. While Lattice QCD (LQCD) calculations of simple hadronic systems are now 
being performed at the physical light-quark masses and the effects of quantum electrodynamics 
(QED) are beginning to be included (see, e.g. Ref. [1]), such calculations have not yet been pre¬ 
sented for more complex systems such as nuclei. However, remarkable progress has been made in 
the ongoing efforts to calculate the lowest-lying energy levels of the simplest nuclei and hypernuclei 
(with A < 4) and the nucleon-nucleon scattering S-matrix elements [2H22] • The magnetic moments 
and polarizabilities of the light nuclei have recently been calculated 123 EH, and by determining 
the short-range interaction between nucleons and the electromagnetic field, the first LQCD calcu¬ 
lation of the radiative capture process np dy [25] was recently performed and the experimentally 
measured cross section was recovered within the uncertainties of the calculation after extrapolation 
to the physical quark masses. These calculations represent crucial steps toward verifying LQCD as 
a useful technique with which to calculate the properties of nuclear systems. However, it will take 
significant computational resources to reduce the associated uncertainties below those of experi¬ 
ment. Near term advances in the field will come from calculations of quantities that are challenging 
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or impossible to access experimentally, such as multi-nucleon forces, hyperon-nucleon interactions, 
rare weak matrix elements and exotic nuclei, such as hypernuclei and charmed nuclei, that are of 
modest computational complexity. Further, performing calculations specifically to match LQCD 
results to low-energy effective field theories (EFTs) will provide a means to make first predictions 
at the physical quark masses and to make predictions of quantities beyond those calculated with 
LQCD. Such calculations are now underway, using the results of our previous works and those 
of Yamazaki et al.^ with the first efforts described, for example, in Ref. |26| for hyperon-nucleon 
interactions and Ref. [271EU for nucleon-nucleon interactions and light nuclei. 

In this work, we present the results of LQCD calculations of two-nucleon systems performed at 
a pion mass of 171^^ ^ 450 MeV in three lattice volumes of spatial extent L — 2.8 fm, 3.7 fm and 
5.6 fm at a lattice spacing of 5 ^ 0.12 fm. In Section [11} we introduce the LQCD methods that are 
used to determine correlation functions and Section |III| reports the results of precision studies of 
the single hadron systems. Section [fv| explores the coupled-channel systems in detail, while 

the channel is discussed in Section]^ In Section |V^ these channels are further investigated in 
the context of nucleon-nucleon effective field theory (NNEFT) before the conclusions of the study 


are presented in Section VII 


II. METHODOLOGY 
A. Calculational Details 

LQCD calculations were performed on three ensembles of nj = 2 -h 1 isotropic gauge-field con¬ 
figurations with L = 24, 32 and 48 lattice sites in each spatial direction, T = 64, 96, 96 sites in the 
temporal direction, respectively, and with a lattice spacing of 6 = 0.1167(16) fm [29]. The Liischer- 
Weisz gauge action 1301 was used with a clover-improved quark action m with one level of stout 
smearing (p = 0.125) [32]. The clover coefficient was set equal to its tree-level tadpole-improved 
value, a value that is consistent with an independent numerical study of the nonperturbative csw 
in the Schrodinger functional scheme [33ll35] . reducing discretization errors from 0{b) to 
The L = 24,32 and 48 ensembles consist of 3.4 x 10^, 2.2 x 10"^, and 1.5 x 10^ HMC evolution 
trajectories, respectively. Calculations were performed on gauge-field configurations taken at uni¬ 
form intervals from these trajectories, see Table [I| The strange-quark mass was tuned to that 
of the physical strange quark, while the selected light-quark mass gave rise to a pion of mass 
rrij^ = 449.9(0.3) (0.3)(4.6) MeV and a kaon of mass mx = 595.9(0.2)(0.2)(6.1) MeV. Many details 
of the current study mirror those of our previous work at the SU(3) symmetric point, which can be 
found in Refs. iniiiT]. In each run on a given configuration, 48 quark propagators were generated 
from uniformly distributed Gaussian-smeared sources on a cubic grid with an origin randomly se¬ 
lected within the volume. Multiple runs were performed to increase statistical precision and the 
total number of measurements is recorded in Table [I] Specifics of the ensembles and the number of 
sources used in each ensemble can also be found in Table |T] Quark propagators were computed using 
the multigrid algorithm [36] or using GPUs jSZlEH] with a tolerance of 10 in double precision. 
In the measurements performed on the L — 24 and 32 ensembles, the quark propagators, either 
unsmeared or smeared at the sink using the same parameters as used at the source, provided two 
sets of correlation functions for each combination of source and sink interpolating fields, labeled 
as SP and SS, respectively. In contrast, for the measurements performed on the L = 48 ensemble 
only SP correlation functions were produced. The propagators were contracted into baryon blocks 
that were projected to a well-defined momentum at the sink, that were then used to form the one- 
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TABLE I: Parameters of the ensembles of gauge-field configurations and of the measurements used in this 
work. The lattices have dimension x T, a lattice spacing 6, and a bare quark mass b rriq (in lattice units). 
Nsrc light-quark sources are used (as described in the text) to perform measurements on N^g configurations 
in each ensemble. 


Label 

A 

B 

C 


L/b T/b /3 b mi b ms b [fm] L [fm] T [fm] m^^L m^^T A'cfg N^rc 
24 64 6.1 -0.2800 -0.2450 0.1167(16) 2.801(29) 7.469(77) 6.390 17.04 4407 1.16 x 10^ 

32 96 6.1 -0.2800 -0.2450 0.1167(16) 3.734(38) 11.20(12) 8.514 25.54 4142 3.95 x 10^ 

48 96 6.1 -0.2800 -0.2450 0.1167(16) 5.602(58) 11.20(12) 12.78 25.49 1047 6.8 x 10^ 


and two-nucleon correlation functions. The blocks are of the form 








(/3)A'/ 


l{N) 


( 1 ) 


where is a quark propagator of flavor f — and the indices are combined spin-color indices 
running over i = 1,..., where A^c = 3 is the number of colors and A^^ = 4 is the number of 

spin components. The choice of the fi and the tensor b^^^ depend on the spin and flavor of the 
nucleon under consideration, and the local interpolating fields constructed in Ref. [39], restricted 
to those that contain only upper-spin components (in the Dirac spinor basis) are used. This choice 
results in the simplest interpolating fields that also have good overlap with the nucleon ground 
states (from localized sources). Blocks are constructed for all lattice momenta |p|^ < 5 allowing 
for the study of two-nucleon systems with zero or nonzero total momentum. In the production on 
the L = 32 ensemble, correlation functions were produced for all of the spin states associated with 
each nuclear species. However, only one spin state per species was calculated on the L — 2A and 
L — A8 ensembles. 


B. Robust Estimators: The Mean with Jacknife and the Hodges-Lehmann Estimator with 

Bootstrap 

The correlation functions are estimated from calculations performed from many source locations 
on many gauge-field configurations. On any given configuration, these results are correlated and, 
because they become translationally invariant after averaging, they can be blocked together to 
generate one representative correlation function for each configuration. More generally, because of 
the correlation between nearby configurations produced in a Markov chain, the results obtained over 
multiple gauge-field configurations are blocked together to produce one representative correlation 
function from any particular subsequence of the Markov chain. In this work, there are a large 
number of independent representative correlation functions which, by the central limit theorem, 
tend to possess a Gaussian-distributed mean. As computational resources are finite, only a finite 
number of calculations of each correlation function can be performed. The underlying distributions 
of the nuclear correlation functions are non-Gaussian with extended tails, and therefore outliers 
are typically present in any sample which lead to slow convergence of the mean. This can then lead 
to significant fluctuations in estimates of correlation functions when resampling methods, such as 
Bootstrap and Jacknife, are employed using the mean to estimate average values (for a discussion 
of the “noise” associated with these and other such calculations, see Ref. pOIIi^ ). Dealing with 
outliers of distributions is required in many areas beyond LQCD, and there is extensive literature 
on robust estimators that are resilient to the presence of outliers, such as the median or the 
Hodges-Lehmann (HL) estimator [43] . The vacuum expectation values of interest in quantum field 
theory are defined by the mean value of a (generally non-Gaussian) distribution. Nevertheless, 
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with sufficient blocking, the mean of the distribution will be Gaussian distributed, for which the 
mean, median, mode and HL estimator coincide. It therefore makes sense to consider such robust 
estimators for large sets of blocked LQCD correlation functions. 

While the median of a sample {xi} is well known, the HL-estimator is less so. It is a robust 
and unbiased estimator of the median of a sample, and is defined as [13] 

HL({xJ) = Median [{(x^ + Xj)/2}] , (2) 

where the sample is summed over all 1 < < A/^, where N is the sample size. The uncertainty 

associated with the HL-estimator is derived from the Median Absolute Deviation (MAD), defined 
as 


MAD({x^}) = Median [{I— Median [{x^}] I}] . (3) 

For a Gaussian distributed sample, la = 1.4826 MADs. The median, HL-estimator and other sim¬ 
ilar estimators cannot be computed straightforwardly under Jackknife, and instead such analyses 
are performed with Bootstrap resampling. 

In the present work, the correlation functions, and their ratios, are analyzed using both the 
mean under Jackknife and HL under Bootstrap, from ^ 100 representative correlation functions 
constructed by blocking the full set of correlation functions. In almost all cases, the HL with 
Bootstrap gives rise to smaller statistical fluctuations over the resampled ensembles and, conse¬ 
quently, to smaller uncertainties in estimates of energies, as seen in our previous investigation into 
robust estimators mi. It is found that outlying blocked correlation functions cause a significant 
enlargement of the estimated variance of the mean, while the robust HL-estimator is insensitive to 
them. 


III. SINGLE MESONS AND BARYONS 

Precision measurements of the single hadron masses, their dispersion relations and their volume 
dependence are essential for a complete analysis of multi-nucleon systems, in particular for a 
complete quantification of the uncertainties in binding energies and S-matrix elements. Single 
hadron correlation functions for the tt^, the octet baryons and the decuplet baryons 

were calculated in each of the three lattice volumes at six different momenta (in each volume), 
from which ground-state energies for each momentum were extracted. The hadron energies were 
extracted from plateaus in the effective mass plots (EMPs) derived from linear combinations (in 
the L = 24 and 32 ensembles) of the SP and SS correlation functions calculated at each lattice 
momentum. The EMPs associated with the tt^ and are shown in Fig.[l]and Fig.[^ respectively, 
while the EMPs for the octet baryons are shown in Figs. iii an dH It is clear from the EMPs 



FIG. 1: Cosh EMPs for the tt^ in the L = 24 (left), L = 32 (center), L = 48 (right) lattice volumes, 
respectively. In ascending order, the momenta are P = 27rn/I/ with |np = 0,1,2, 3,4, 5. 
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FIG. 2: Cosh EMPs for the in the L = 24 (left), L = 32 (center), L = 48 (right) lattice volumes, 
respectively. In ascending order, the momenta are P = 27rn/I/ with |np = 0,1,2, 3,4, 5. 




FIG. 3: EMPs for the nucleon in the L = 24 (left), L = 32 (center), L = 48 (right) lattice volumes, 
respectively. In ascending order, the momenta are P = 27rn/I/ with |np = 0,1,2, 3,4, 5. 




FIG. 4: EMPs for the A in the L = 24 (left), L = 32 (center), L = 48 (right) lattice volumes, respectively. 
In ascending order, the momenta are P = 27rn/I/ with |np = 0,1, 2, 3,4, 5. 




FIG. 5: EMPs for the E in the L = 24 (left), L = 32 (center), L = 48 (right) lattice volumes, respectively. 
In ascending order, the momenta are P = 27rn/I/ with |np = 0,1, 2, 3,4, 5. 


that extended ground-state plateaus exists for all hadrons at all momenta, and as such relatively 
precise hadron masses and dispersion relations can be determined. A correlated x^-i^iiii^ization 
fit of the plateau region in combinations of correlation function to a constant energy was performed 
over a range of fit intervals to determine the energy, its statistical uncertainty and the systematic 
uncertainty due to the selection of the fitting range. The energies of the pseudoscalar mesons and 
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FIG. 6: EMPs for the E in the L = 24 (left), L = 32 (center), L = 48 (right) lattice volumes, respectively. 
In ascending order, the momenta are P = 27rn/I/ with |np = 0,1, 2, 3,4, 5. 


octet baryons are are shown in Table [IT] and Table in respectively. 


TABLE II: The pion and kaon energies (l.u.) as a function of momentum (l.u.), |P| = (^) |n|, calculated 
on each ensemble of gauge-field configurations. The infinite-volume meson masses, determined by fitting 
expressions of the form in eq. are also given. The first uncertainty associated with each extraction is 
statistical and the second is the fitting systematic. In the case of the extrapolated values, the systematic 
uncertainty also contains the estimated uncertainty due to the extrapolation (which is small in both cases). 


meson 

ensemble 

|n|=0 |n|2 = l |n|2 = 2 |n|2 = 3 |n|2 = 4 |n|2 = 5 

TT^ 

243 X 64 
32^ X 96 
48^ X 96 

0.26626(36)(14) 0.37184(28)(34) 0.45341 (30)(45) 0.5204(07)(13) 0.5812(10)(17) 0.6329(09)(12) 
0.26607(23)(09) 0.33006(20)(14) 0.38330(21)(16) 0.43042(26)(28) 0.47156(43)(93) 0.5093(05)(12) 
0.26607(17)(11) 0.29624(14)(05) 0.32365(13)(10) 0.34895(16)(10) 0.37221(22)(18) 0.39404(31)(35) 


L = 00 

0.26606(14) (08) 

K± 

243 X 64 

323 ^ gg 

483 X 96 

0.35239(30)(16) 0.43749(24)(25) 0.50810(22)(25) 0.56947(35)(50) 0.6224(07)(13) 0.67109(52)(55) 
0.35248(18)(08) 0.40259(16)(17) 0.44725(17)(09) 0.48782(24)(49) 0.52357(45)(60) 0.55727(46)(88) 
0.35236(16)(25) 0.37559(13)(06) 0.39744(13)(06) 0.41814(13)(06) 0.43760(17)(05) 0.45628(21)(09) 


L = 00 

0.35240(11)(03) 


The energies determined at zero momentum are used to extrapolate the hadron masses to 
infinite volume, and are combined with the other energies to determine their dispersion relations. 
With the large values of m^L in the ensembles of gauge configurations, it is sufficient to use the 
leading-order (LO) finite-volume (FV) corrections to the hadron masses to extrapolate from the 
volumes of the calculations to infinite volume. The LO modifications to the pseudoscalar masses, 
tum^ and baryon masses, Mb^ are given by 

iy)( T\ M , 

[m^L) = m>M + - 

M^Pim^L) = ^ + ... , (4) 

where the forms are those of p-regime chiral perturbation theory (xPT) and heavy-baryon xPT 
(HBxPT [45]). The infinite-volume masses, and M^\ and the coefficients of the LO volume 
dependence, cm and c^, are quantities determined by fits to the LQCD calculations, and will, in 
general, be different for each hadron. 

The zero-momentum energies of the pseudoscalar mesons and their infinite-volume extrapolation 
are given in Table |TT| and shown in Fig. The energies of both mesons are found to be independent 
of the lattice volume within the uncertainties of the calculations. Despite the larger number of 
correlation functions in the L = 24 ensemble, the uncertainties in the meson masses are larger than 
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TABLE III: The baryon energies (l.u.) as a function of momentum (l.u.), |P| = (^) |n|, calculated on each 
ensemble of gauge-field configurations. The infinite-volume masses, determined by fitting the expression in 
eq. are also given. The first uncertainty is statistical and the second is the fitting systematic. In the 
case of the extrapolated values, the systematic uncertainty also contains the estimated uncertainty due to 
the extrapolation (which is small in all cases). 


baryon 

ensemble 

II 

o 

to 

II 

to 

II 

to 

to 

II 

CO 

to 

II 

to 

II 

or 

N 

24® X 64 
32=* X 96 
48=* X 96 

0.7251(04)(11) 0.7699(10(13) 0.8108(10)(13) 0.8497(13)(21) 0.8944(16)(23) 0.9311(17)(23) 

0.72546(47)(31) 0.75160(60)(47) 0.77657(75)(89) 0.80098(62)(81) 0.8238(07)(11) 0.8467(07)(10) 
0.7245(10)(13) 0.7359(21)(34) 0.7471(23)(35) 0.7556(20) (36) 0.7661(21)(40) 0.7771(22)(42) 


L = oo 

0.72524(46) (35) 

A 

24=* X 64 
32=* X 96 
48=* X 96 

0.77609(42)(66) 0.8165(14)(18) 0.8533(14)(21) 0.8918(23)(34) 0.9336(14)(22) 0.9709(12)(16) 

0.77633(45)(48) 0.80059(60)(48) 0.82435(75)(51) 0.84687(78)(54) 0.8680(10(14) 0.8900(08)(10) 

0.77650(94)(80) 0.7858(14)(20) 0.7963(14)(21) 0.8066(15)(23) 0.8166(16)(27) 0.8268(16)(29) 


L = OO 

0.77638(42)(48) 

E 

24=* X 64 
32=* X 96 
48=* X 96 

0.79520(70)(65) 0.83608(73)(62) 0.87550(75)(87) 0.9147(07)(13) 0.9485(07)(10) 0.9855(10)(24) 
0.79634(31)(49) 0.82033(60)(61) 0.84320(63)(75) 0.86502(71)(51) 0.88575(60)(57) 0.90755(65)(63) 
0.7958(12)(13) 0.8050(14)(23) 0.8152(15)(24) 0.8253(16)(26) 0.8351(16)(28) 0.8451(17)(30) 


L = OO 

0.79638(33)(54) 


24=* X 64 
32=* X 96 
48=* X 96 

0.83646(63)(49) 0.87594(60)(58) 0.91318(58)(54) 0.9487(06)(10) 0.9828(06)(11) 1.01668(60)(95) 
0.83715(53)(58) 0.85886(49)(59) 0.88044(50)(57) 0.90201(51)(36) 0.92261(62)(89) 0.94276(66)(89) 
0.83643(68)(72) 0.8460(11)(10) 0.8557(12)(11) 0.8652(12)(13) 0.8744(13)(14) 0.8837(14)(17) 


L = OO 

0.83690(45)(50) 





FIG. 7: The volume dependence of the 7r+ (left panel) and (right panel) masses. Energies (l.u.) in 
the L = 24, 32 and 48 lattice volumes are shown as the blue, yellow and red points, respectively, while the 
results of fits to these results of the form given in eq. Q are shown by the shaded regions with the inner 
(outer) band denoting the statistical (statistical and systematic combined in quadrature) uncertainties. 


those extracted from the L = 32 ensemble. The zero-momentum energies of the octet baryons 
and their infinite-volume extrapolation are given in Table m and shown in Fig. As with the 
mesons, there is no statistically significant volume dependence observed for any of the octet-baryon 
masses. Two-parameter x^-i^ii^imization fits of the form given in eq. Q were performed to the 
volume-dependence of each hadron to extract its infinite-volume mass. Because of the negligible 
volume dependence in the LQCD results, limited constraints can be placed on the cm,b coefficients. 
In addition to the tt^, and octet baryons, analogous extrapolations were performed with the 
results obtained for the and decuplet baryons, the combined results of which are shown 
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FIG. 8: The volume dependence of the A, H and E masses. Energies (l.u.) in the L = 24,32 and 
48 lattice volumes are shown as the blue, yellow and red points, respectively, while fits to these results 
are shown by the gray, shaded regions with the inner (outer) band denoting the statistical (statistical and 
systematic combined in quadrature) uncertainties. 


in Table [Tv| (in both l.u. and MeV). 


TABLE IV: The infinite-volume hadron masses obtained by extrapolating zero-momentum ground-state 
energies with the volume dependence given in eq. ©• The first and second uncertainties are the statistical 
and systematic, respectively, while the third for values in units of MeV results from the uncertainty in the 
scale setting. 


hadron 

Mass (l.u.) Mass (MeV) 

hadron 

Mass (l.u.) Mass (MeV) 

TT^ 

0.26614(15)(15) 449.9(0.3)(0.3)(4.6) 


0.35241(12)(11) 595.9(0.2)(0.2)(6.1) 


0.5248(14)(15) 887.3(2.4)(2.5)(9.1) 

K*± 

0.56923(89)(51) 962.4(1.5)(0.9)(9.9) 

N 

0.72524(46) (35) 1226(01)(01)(12) 

A 

0.77638(42) (48) 1312(01)(01)(13) 

E 

0.79638(33) (54) 1346(01)(01)(14) 


0.83690(45)(50) 1415(01)(01)(15) 

A 

0.8791(14)(17) 1486(02)(03)(15) 

E* 

0.9211(17)(19) 1557(03)(03)(16) 


0.9637(09) (17) 1629(02)(03)(17) 


1.0059(06)(12) 1700(01)(02)(17) 


Deviations of the single hadron dispersion relations from that of special relativity lead to mod¬ 
ifications to Liischer’s quantization conditions (QCs) in two-body systems. To address this, the 
dispersion relations have been precisely determined, and the deviations from special relativity are 
propagated through the extraction of S-matrix elements using the QCs. In each of the ensembles, 
single hadron correlation functions were calculated for each of the hadrons of interest with mo- 
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menta |p| < \/5 jL\ the results of which are given in Table [n| and Table III Energy-momentum 
relations that are fit to the results obtained for each hadron, /i, are of the form 




(5) 


where the hadron speed of light, and the higher-order deviation from special relativity, pa¬ 
rameterized by %, are determined by fits to the results of the LQCD calculations. With this 
parameterization, the V}^ are consistent with unity and the r]h sire consistent with zero (for all 
hadrons). There is a Lorentz-breaking term that could be considered at this order in a momentum 
expansion, but this is also found to be consistent with zero. The energies of the tt^, 

3 




FIG. 9: Dispersion relations of the tt^, . The results in the L = 24, 32 and 48 lattice volumes are shown 
as the blue, yellow and red points, respectively, while fits to these results are shown by the gray curves. 


and octet baryons as a function of momentum are given in Table [IT] and Table III and shown in 
Fig. and Fig. W X^-minimization fits to the energy-momentum dispersion relation are performed 


to extract the speed of light for each hadron, the results of which are shown in Table jVj In the 


TABLE V: The speed of light of each hadron determined from fits to the energy-momentum results. 


hadron 

Vh 

hadron 

Vh 

TT^ 

1.0025(18)(08) 


1.0038(20)(12) 

N 

1.010(16)(07) 

A 

1.018(15)(01) 

E 

1.010(12)(03) 

5 

1.0102(61)(13) 


low-energy regime relevant to the two-nucleon systems, the dispersion relation of special relativity 
is found to hold at the ^ 1% level. 


A. The Gell-Mann-Okubo Mass Relation 

Given the precise determinations of the single hadron spectrum, it is important to test relations 
between baryon masses that are predicted to hold in particular limits of QCD. The Gell-Mann- 
Okubo mass relation [SI ST] arises from SU(3) flavor symmetry and its violation, quantifled by 

Tqmo — 


( 6 ) 
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(27Tb/L)^n2 (27Tb/L)^n2 




(277b/L)^ (277b/L)^ ri^ 


FIG. 10: Dispersion relations of the octet baryons. The results in the L = 24, 32 and 48 lattice volumes 
are shown as the blue, yellow and red points, respectively, while fits to these results are shown by the gray 
curves. 


results from SU(3) breaking transforming in the 27-plet irreducible representation (irrep) of flavor 
SU(3) which can only arise from multiple insertions of the light-quark mass matrix or from non- 
analytic meson-mass dependence induced by loops in xPT. Further, it has been shown that Tgmo 
vanishes in the large-A^c limit as 1/Nc [48]. In previous work [49], we performed the first LQCD 
determination of this quantity, after which more precise LQCD determinations m were performed. 
In this work, by far the most precise determination of Tqmo was obtained from the L = 32 
ensemble, where we find Tgmo = +0.000546(51)(81) l.u. = -h0.92(09)(14)(01) MeV (compared with 
Tgmo = +0.00056(19)(38) l.u. = +0.96(33)(64)(01) MeV and Tgmo = +0.00104(27)(29) l.u. = 
+1.76(46)(49)(02) MeV, from the L = 24 and 48 ensembles, respectively). It is conventional to 
form the dimensionless quantity ^gmo = Tqmo/Mq^ where Mq is the centroid of the octet baryons 
masses. In the present calculations, the centroid is found to be Mq = 0.78658(51) (36) l.u. = 
1329(01)(01)(14) MeV, from which ^gmo = 0.00069(06)(10). This value is consistent with our 
previously published result close to this pion mass and is also consistent with other subsequent 
determinations [50], but far more precise. It is worth noting that the experimental value, = 

+8.76(08) MeV, is an order of magnitude larger than the value we have determined at this heavier 
pion mass. 


IV. THE COUPLED CHANNELS AND THE DEUTERON 

The phenomenology of the ^Si-^Di coupled J = 1 channels in finite volumes has been explored 
recently using the experimentally constrained phase shifts and mixing angles in an effort to under¬ 
stand what might be expected in future LQCD calculations m- One goal of that study was to 
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estimate the lattice volumes, and identify the correlation functions, required to extract the phase 
shifts and mixing parameter describing these channels in infinite volume. It was found to be conve¬ 
nient in those FV studies m to use the Blatt-Biedenharn (BB) [52] parameterization of the 2x2 
S-matrix (below the inelastic threshold). 




( cos Cl — sin ei 
sin ei cos ei 



( cos 6i sin 6i 
— sin ei cos ei 


(7) 


from which the QCs associated with these channels can be determined. For the two-nucleon system 
at rest in a cubic volume, embedded in the even parity Ti irrep of the cubic group, the QC in the 
limit of vanishing 5iy3, D-waves and higher phase shifts becomes [5T] 

= 47rcS’°’°)(fc|^;L) , (8) 


where is the magnitude of the momentum in the center-of-momentum (CoM) frame, and the 

function ; L) is proportional to the Liischer Zqq function, as given in Ref. [53l[5l]. The 

phase shift 6ia is evaluated at . The three j^-substates are degenerate and their energies are 
insensitive to the mixing parameter 6i. 

In contrast, for the two-nucleon system carrying one unit of lattice momentum along the z-axis, 
Ptot. = with d = (0,0,1), the three substates are embedded into two distinct even-parity 
irreps of the cubic group - the one-dimensional A 2 representation and the two-dimensional E 
representation, containing the jz = 0 and jz = ±1 states, respectively. In the same limit as taken 
to derive eq. the QCs for these two irreps are m 


fcla cot Sia{kl^) 


k^COtSi^ik^) 


47rcS^^'^\kl^] L) - L) Se.ikl^) 


(9) 


where 


= v^sin26i(A:*) — sin^ 6 i(A;*) . (10) 

The difference in energy between the A 2 and E FV eigenstates provides a measure of ei, but this 
is complicated by the fact that they are evaluated at two slightly different energies. This analysis 
can be extended to other lattice momenta m, but the QCs in eq. Q and eq. are sufficient for 
the present purposes. 

Correlation functions for two nucleons in the Ti, A 2 and E irreps are straightforwardly con¬ 
structed from the nucleon blocks we have described previously. In fact, multiple correlation func¬ 
tions are constructed in each irrep. In the L — 24 and 48 ensembles, the spin projections were not 
performed to permit construction of the A 2 irrep, and so only L = 32 correlation functions can be 
used to constrain ei. 


A. The Deuteron 

In nature, the deuteron is the only bound state in the two-nucleon systems, residing in the 
coupled channels, and it has a special position in nuclear physics. The deuteron has always 
provided a benchmmark when deriving phenomenological interactions between nucleons, and it 
will play a critical role in verifying LQCD as a viable calculational tool. Correlation functions for 
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TABLE VI: The deuteron binding energies extracted from plateaus in the EMPs shown in Eig. [TO along 
with the infinite-volume extrapolated value. The size of the EV effects is characterized by shown 

in the last column. The first uncertainty corresponds to the statistical uncertainty associated with the 
fit, the second corresponds to the systematic uncertainty associated with the selection of the fitting interval 
(determined by varying this range). In the case of dimensionful quantities, the third uncertainty is associated 
with scale setting. Eor the infinite-volume values of the binding energy, the last uncertainty is introduced 
by the finite-volume extrapolation in eq. 0, and is estimated by considering the effect of omitted terms 
scaling as jL. 


Ensemble 

AE (l.u.) 

Bd (MeV) 

g-«L 

243 X 64 

-0.01157(73) (96) 

19.6(1.2)(1.6)(0.2) 

0.111 

32^ X 96 

-0.01037(89) (96) 

17.5(1.5)(1.6)(0.2) 

0.063 

48^ X 96 

-0.0078(12)(19) 

13.3(2.0)(3.2)(0.2) 

0.027 

L = 00 


, . .-K(1.6)(2.7)(0.2)(0.2) 
-^^•^-(1.8)(1.8)(0.2)(0.2) 



two nucleons in the even-parity Ti irrep of the cubic group were constructed, from which, after 
a correlated subtraction of twice the energy of a single nucleon, the EMPs shown in Fig. [TT] were 
derived As with the single hadrons, correlated x^-i^iinimization fits of a constant to the plateau 


0.01 

0 


- - 0.01 
- 0.02 

0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 

t/b t(l.u.) t(l.u.) 

EIG. 11: EMPs for the energy difference between the deuteron and twice the nucleon in the L = 24 (left), 
L = 32 (center) and L = 48 (right) ensembles, along with fits to the plateau regions. The extracted binding 
energies are given in Table |VT[ 












i 



regions were performed to estimate the deuteron binding energy, and associated uncertainties. 
The deuteron binding energies extracted from each ensemble are given in Table VI, along with the 
values of where n — yjMjsiBd is the binding momentum of the deuteron. As is seen to 

change from ^ 3% in the largest volume to ~ 11% in the smallest, an extrapolation in volume is 
desirable. 

Inspired by the FV contributions to the binding of a shallow bound state resulting from short- 
range interactions [SSHSZ], the extrapolation to infinite volume was performed by fitting a function 


^ The single nucleon correlation function, the square of which is divided out of two-nucleon correlation functions 
to yield a plateau on the energy difference, have been temporally displaced, in some instances, to enhance the 
plateau region in the difference. Further, due to the nature of the HL estimator, the first few time-slices in the 
difference correlation function have been removed, leading to temporal displacements of the EMPs. The EMPs 
defining energy differences in this work correspond to both the one-nucleon and two-nucleon correlation functions 
being in their respective ground states (as defined by plateaus in their respective individual EMPs). 
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of the form, 


Bd{L) = + Cl 


^-\/2koL ^ ^-VSkqL 

—— + V2 --- + ^--- 

L L 3\/3 L 


+ 


( 11 ) 


to the results obtained in the three lattice volumes, where kq = \/ (with B'^’ the deuteron 

binding energy in infinite volume) and ci are the fit parameters. The ellipsis denote terms that 
are and higher. A x^-minimization fit to the deuteron binding energies in Table VI 

" +1. The 


?(oo) 


generates the region in parameter space shown in Fig. 12, defined by Xmin 



C1 


FIG. 12: The region in parameter space defined by ^ Xmin + inner region is defined by 

the statistical uncertainty, while the outer region is defined by the statistical and systematic uncertainties 
combined in quadrature. 


deuteron binding energy found from extrapolating to infinite volume is 


- 14 4 +( 1 . 6 )( 2 . 7 )( 0 . 2 )( 0 . 2 ) 


( 12 ) 


The first uncertainty corresponds to the statistical uncertainty, the second corresponds to the 
fitting systematic uncertainty, the third is associated with scale setting, and the last uncertainty 
is introduced by the finite-volume extrapolation in eq. and is estimated by considering the 

effect of terms scaling as ^ Combining the errors in eq. (12) in quadrature leads to 

= UAtli MeV. 


1. The Mixing Parameter, ei 


For a deuteron that is moving in the lattice volume, the energy eigenvalues are sensitive to the 
mixing parameter ei, as expected from the QCs given in eq. Q for the specific boost d = (0,0,1). 
Explicitly evaluating the cf^ functions that appear in eq. for the two irreps containing the 
deuteron gives the QCs, 


^A2 ^1q;('^^A2) T ^A2 

cot 5la(mE) + = 


2e 


— HA2 ^ 


L 

2e-«Ei 


1 + 2 1 + 


+ 


KA2L (kajI/)^ 


Sji {inpi^2 


1-1 + 


+ 


hkL (ke-C')^ 




(13) 


where Se-^{k*) is defined in eq. (10). For both irreps, the functions k*cotS and are evaluated 
at k* = in. Iteratively solving these QCs in terms of the infinite-volume binding momentum, kq 
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(ka 2 ) in the infinite-volume limit), the spin-averaged binding energy of the A 2 and E irreps 

is 


r(0,0.1) 


= B 


(00) 


-«oL , 
ML ^ 


(14) 


where the ellipses denote terms and higher, which is consistent, at this order, with the 

binding energy extracted from the Ti irrep for the deuteron at rest. In the above expression, 
is the residue of the deuteron pole. The difference in energies is 


5B 


( 0 , 0 , 1 ) _ 


12koZI 




M L 


1 + 


+ 


KqL 


^ei (^^ 0 ) “ 1 “ 


(15) 


Calculating the exponentially small difference between the energies of these two states provides a 
direct measure of ei evaluated at the deuteron pole. In order to extract a meaningful constraint 
on Cl, the FV corrections must be statistically different from zero, otherwise the coefficient of the 
leading contribution to the energy difference vanishes. 

In the present production, it has been only possible to decompose the d = (0,0,1) boosted 
deuteron correlation functions into the E (j^ = ±1) and A 2 {jz — 0) irreps in calculations performed 
with the L — 32 ensemble. The EMP associated with the difference in energies between these irreps 


0.03 

_ 0.020 
q 0.010 


< - 0.010 
- 0.020 
- 0.03 

0 2 4 6 8 10 12 14 

t(l.u.) 

FIG. 13: EMP associated with the energy difference between the E {jz = ±1) and A 2 {jz = 0) deuteron 
states with boost vector d = (0,0,1) in the L = 32 ensemble, along with fits to the plateau region. The 
energy difference depends upon the mixing parameter ei. 



is shown in Fig. and the energy difference extracted from fitting the plateau region is consistent 
with zero, {L = 32) = —0.4(4.1) (4.6) MeV. While this energy difference is bounded in 

magnitude, the fact that the FV contributions to the deuteron binding energy are consistent with 
zero in this lattice volume means that no useful bound can be placed upon ei. 


2. A Compilation of Deuteron Binding Energies from LQCD 


The current calculation of the deuteron binding energy adds to a small number of previous 
calculations over a range of pion masses above ^ 300 MeV [HI [HI [151 [2T] , ^ shown in Fig. 


14 


^ The deuteron and dineutron binding energies at tUtt ~ 800 MeV in the L = 24 and L = 32 ensembles presented 
in Ref. [TJ have been reproduced in Ref. |22j . within uncertainties, on the same gauge ensembles. 

^ The results of quenched calculations, and of calculations that have not been extrapolated to infinite volume [2], have 
not been shown. The results from Ref. [151121] were obtained with a power-law extrapolation to infinite volume. 
This is not the correct form for a loosely bound state, and tends to lead to signihcantly smaller uncertainties than 
from extrapolations performed with the known exponential form. 
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The present result is consistent, within uncertainties, with the results at ^ 300 MeV and 
^ 500 MeV from Refs. [T5l |2T]. Further LQCD calculations at lighter quark masses are 


QQ 


30 

25 

.20 

10 


) NPLQCD, anisotropic 
) Yamazaki et al. 

) NPLQCD, isotropic 


♦ ( 


200 400 600 

rrijj (MeV) 


800 


FIG. 14: The pion mass dependence of the denteron binding energy calculated with LQCD. The NPLQCD 
anisotropic-clover result is from Ref. m , the Yamazaki et al results are from Refs. [151 [21] and the NPLQCD 
isotropic-clover results are from this work and Ref. m- The black disk corresponds to the experimental 
binding energy. 


required to quantify the approach to the physical denteron binding (for related NNEFT work see 
Ref. [28]b 


B. Scattering in the Coupled Channels 


To recover the S-matrix in the coupled channels, calculations must be performed that 

isolate the phase shifts and mixing angle, 5ic^, ei and defined in eq. Q, from the FV observables 
accessible to LQCD calculations. The formalism with which to perform this analysis [5ll|58H6n| 
is an extension of the seminal work of Liischer [531 154j. For vanishing total momentum, assuming 
that the contribution from D-waves and higher are negligible, the energies of the Ti irreps are 

insensitive to ei, as demonstrated in eq. Q. Therefore, the shifts in energies of the two nucleon 
states in the Ti irrep for various total momentum from the energy of two free nucleons can be used 
to extract 6ia below the inelastic threshold. 

Figure 15 show the effective-/c*^ plots (Ek2Ps) associated with the first continuum Ti states 
in each ensemble, with momentum near k = 271 jL. These show the values of the interaction 
momentum extracted from the LQCD correlation functions as a function of Euclidean time. 
As with the EMPs, plateau behavior indicates the dominance of a single state. For an arbitrary 



t(l.u.) t(l.u.) t(l.u.) 


FIG. 15: Ek2Ps for the lowest lying continuum NN states near /c* = 271 /L in the L = 24 (left), 

L = 32 (center) and L = 48 (right) ensembles, along with fits to the plateau regions. 


two-body system, comprised of particles with masses mi and m 2 , with zero CoM momentum, the 
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interaction momentum is defined through 


= S’* — mi — m2 = ^ A:*^ + m^ — mi — m2 


(16) 


where £’* is the energy in the CoM frame, defined by £’* = — |Ptot.p where E is the total 

energy, and Ptot. is the total momentum, of the system. Figure [T^ shows the Ek2Ps for states 


with momentum near k = 4:7tIL^ while Fig. 17 shows the Ek2P for the system with d = (0,0,1) on 
the L = 32 ensemble. Inserting the values of fc* extracted from the plateau regions of the Ek2Ps 



t(l.u.) 


t(l.u.) 


FIG. 16: Ek2Ps for the continuum NN states near k = Air/L in the L = 32 (left) and L = 48 

(right) ensembles, along with fits to the plateau regions. 
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^ 0.01 
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FIG. 17: Ek2Ps for the spin-averaged continuum NN states with d = (0,0,1) near k = 0 in the 

L = 32 ensemble, along with the fit to the plateau region. 


in Fig. and Fig. 16 into the QC in eq. gives rise to the values of A:* cot 6ia and 6ia given 
in Table |Vn| and shown in Fig. Additionally, the result of inserting the value of A:* extracted 
from the plateau in Fig. 17 into the QC for the A 2 and E irreps in eq. is shown in Table |VII| 
and Fig. The uncertainties in each of the extractions are relatively large, magnified by their 


close proximity to a singularity in the kinematic functions Cqq. Even subject to these issues, a zero 
in the phase shift is visible near A:* ~ m^^ ^ 450 MeV, indicative of an attractive interaction with 
a repulsive core. It is interesting to compare this phase shift, at a pion mass of ^ 450 MeV, 
with that of nature, illustrated by the dashed curve in Fig. [T^ The phase shift resulting from a 
partial-wave analysis of experimental data is consistent, within uncertainties, with the phase shift 
calculated at ^ 450 MeV over a large range of momenta. The zeros of the phase shift occur 
at different momenta, but they are nearby. Without results at smaller A:*, a precise extraction of 
the scattering parameters, such as the scattering length and effective range, is not feasible, and 
additional calculations are required in order to accomplish this. However, the determination of 
the binding energy and the two continuum states that lie below the threshold of the t-channel cut 
(set by the pion mass. A:* = m^l2) can be used to perform an approximate determination of the 
inverse scattering length and effective range. A linear fit was performed. A:* cot 6 = — 1/a + 
as shown in Fig. The range of linear fits straddle A:* cot 5 = 0 at A:* = 0, and as such allows 
both = ±oc, and it is useful to consider the constraints on l/a^^*^^) rather than The 
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TABLE VII: Scattering information in the coupled channels. A indicates that the uncertainty 

extends across a singularity of the Liischer function, or that it is associated with the bound state. The 
uncertainties in these quantities are highly correlated, as can be seen from Fig. 


Ensemble 

iPtotl (l.U.) 

k* cot Sia/rn^^ Sia (degrees) 

All 

0 


243 X 64 

0 

0.9754+[^Jg®; - 3.1(1.7)(3.7) 

32^ X 96 

323 X 96 

323 X 96 

0 

0 

I 

0.702+[“;g2 2.3+g;S9) 17(5)(11) 

l-065kS(i?) -11.1(3.8)(8.5) 

„ „„„+(26)(29) „ o,+(24)(15) o„+(13)(23) 

■•■^•'J^-(59)(20) ■'■'^®-(ll)(16) 

483 X 96 

483 X 96 

0 

0 

0.426(03)(12) 0.45+g^;g5 44+gq(08) 

0.662(08)(29) 0.35+;°;JJg;°l) 



'o 



100 200 300 

k (MeV) 


400 


500 


FIG. 18: Scattering in the coupled channels. The left panel shows cot Sia/ra-j^ as a function of 

while the right panel shows the phase shift as a function of momentum in MeV, assuming that 6i^ 
and the D-wave and higher partial-wave phase shifts vanish. The thick (thin) region of each result correspond 
to the statistical uncertainty (statistical and systematic uncertainties combined in quadrature). The black 
circle in the right panel corresponds to the known result from Levinson’s theorem, while the dashed-gray 
curve corresponds to the phase shift extracted from the Nijmegen partial-wave analysis of experimental 
data [6T] . 


correlated constraints on and are shown in Fig. 19 

and effective range determined from the fit region in Fig. [T^ are 


The inverse scattering length 






-1 


-1 


-0.04 


+(0.07)(0.08) 

-(0.10)(0.17) 






„ .+(2.2)(3.5) 
^•»-(1.5)(1.7) 


-0.09 


+(0.15)(0.19) 

-(0.23)(0.39) 


fm ^ 




O .+(1.0)(1.5) 
'^•^-(0.7)(0.8) 


fm 


(17) 


Further calculations in larger volumes (and hence at smaller k*'^) will be required to refine these 
extractions. There is a potential self-consistency issue raised by the size of the effective range that 
is within the uncertainties that are reported. Luscher’s method is valid only for the interaction 
ranges R T/2, otherwise the exponentially small corrections due to deformation of the inter¬ 
hadron forces become large. Assuming the range of the interaction is of similar size to the effective 
range (as expected for ’’natural” interactions), this requirement is not met and deviations from the 
assumed linear fitting function should be entertained. Higher precision analyses will be required 
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FIG. 19: Scattering in the ^Si-^Di coupled channels below the start of the t-channel cut, < m^/4, 
assuming that 6i^ and the D-wave and higher partial-wave phase shifts vanish. The left panel shows 
solid region corresponding to linear fits associated with the statistical uncertainty and the statistical and 
systematic uncertainties combined in quadrature. The right panel shows the scattering parameters, 1/a^ 
and determined from fits to scattering results below the t-channel cut. The solid circle corresponds to 

the experimental values. 


to investigate this further. 


V. THE % CHANNEL AND THE DINEUTRON 

The analysis of LQCD calculations in the ^Sq channel are somewhat simpler than in the ^Si-^Di 
coupled channels as scattering below the inelastic threshold is described by a single phase shift, 
In FV, the relation between energy eigenvalues of the system at rest in the Ai cubic irrep 
and given by eq. ^ with 6ia Unfortunately, the correlation 

functions in this channel have larger fluctuations and excited state contamination than those in 
the ^Si-^Di coupled channels system. Consequently, the uncertainties associated with each energy 
level are larger. 


A. The dineutron 


Unlike in nature, the dineutron is bound at heavier quark masses [H la HU EH EH- Plateaus 
identified with a negatively shifted dineutron were found in all three ensembles, with the associated 
EMPs shown in Fig. and the extracted energy shifts shown in Table VIII[ Performing a volume 



0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 0 2 4 6 8 10 12 14 

t(l.u.) t(l.u.) t(l.u.) 


FIG. 20: EMPs for the dineutron in the L = 24 (left), L = 32 (center) and L = 48 (right) ensembles, along 
with fits to the plateau regions. The extracted binding energies are given in Table |VIII[ 
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TABLE VIII: The dineutron binding energies from fitting to the EMPs shown in Eig. 


Ensemble 

AE (l.u.) 

Bnn (MeV) 

g-KL 

243 X 64 

-0.0142(09)(27) 

24.1(1.5)(4.5) 

0.088 

32^ X 96 

-0.0109(09)(20) 

18.4(1.5)(3.3) 

0.058 

48^ X 96 

-0.0070(11)(18) 

11.8(1.9)(3.1) 

0.033 

L = oc 

^•^^'^-(ll)(27)(01) 

. + (1.7)(2.5)(0.2)(0.2) 
-L^-0(1.9)(4.5)(0.2)(0.2) 



extrapolation using the form given in eq. (11) leads to a binding energy of 

o(oo) _ -12 ^+(1-7)(2.5)(0.2)(0.2) _ ^ 

-°nn — ■^^•'-’-(1.9)(4.5)(0.2)(0.2) 


(18) 


Combining the errors in eq. (18) in quadrature leads to = 12.51H MeV. The ci-B 


doo) 

^nn 



FIG. 21: The region in parameter space defined by Xmin + inner region is defined by 

the statistical uncertainty, while the outer is defined by the statistical and systematic uncertainties combined 
in quadrature. 

parameter space defined by Xmin + ^ determined from an uncorrelated fit to the dineutron 

binding energies in the three volumes is shown in Fig. This dineutron binding energy is 
consistent with the binding energy of the deuteron within uncertainties. The EMPs associated 
with the difference between the deuteron and dineutron energies in each ensemble are shown in 



t(l.u.) t(l.u.) t(l.u.) 


FIG. 22: EMPs for the energy difference between the dineutron and the deuteron in the L = 24 (left), 
L = 32 (center) and L = 48 (right) ensembles, along with fits to the plateau regions. 


^ Extrapolating with a form consistent with a scattering state, which would display a volume dependence of 
AE ~ 1/T^, results in a poor goodness-of-fit. 
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TABLE IX: Energy differences between the dineutron and deuteron from fitting to the EMPs shown in 
Eig. All differences are consistent with zero, as is their infinite-volume extrapolation. 


Ensemble 

Fnn ^deut (hu.) 

E'nn - Edeut (MeV) 

243 X 64 

+0.0022(16)(28) 

+3.7(2.8)(4.7)(0.0) 

323 X gg 

-0.0014(09)(15) 

-2.4(1.6)(2.5) 

483 X 96 

+0.0027(04) (31) 

+4.6(0.7)(5.3) 


extracted. 


1. A Compilation of Dineutron Binding Energies from LQCD 

The current calculation of the dineutron binding energy adds to a small number of previous 
calculations, a compilation of which is shown in Fig. There does not appear to be a clear 

20 
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FIG. 23: The pion-mass dependence of the dineutron binding energy calculated with LQCD. The NPLQCD 
anisotropic-clover result is from Ref. m , the Yamazaki et al. results are from Refs. [T5l[2T] and the NPLQCD 
isotropic-clover results are from this work and Ref. [14]. The black disk corresponds to the location of the 
near-bound state at the physical quark masses. 

pattern emerging as to how the dineutron will unbind as the pion mass is reduced. The results 
that have been obtained in Refs. dSlEI] have consistently smaller uncertainties than those found in 
Ref. [m [13] and in the present work. However, the results are consistent within the uncertainties. 


> NPLQCD, anisotropic 
) Yamazaki et al. 

) NPLQCD, isotropic 


B. Scattering in the ^Sq Channel 


Correlation functions for two nucleons in the ^Sq state were constructed in the Ai irrep of the 
cubic group. The Ek2Ps associated with the states near the k = 27t/L and k = 47r/L noninteracting 

respectively. 


levels are shown in Fig. and Fig. respectively. For the lowest-lying “continuum” state, 
plateaus were found in all three ensembles, however, only the L=32 ensemble has correlation 
functions that were sufficiently clean to extract the next higher level. A plateau was also identified 

The values of k cot and 


26 


in the system with one unit of total momentum, as shown in Fig. 

the phase shift are given in Table and shown in Fig. ^7\ Many of the qualitative features of 


the results for the scattering amplitude in this channel are similar to those in the ^Si-^Di coupled 
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t(l.u.) t(l.u.) t(l.u.) 


FIG. 24: Ek2Ps for the lowest lying continuum ISo NN state near k = 27r/L in the L = 24 (left), L = 32 
(center) and L = 48 (right) ensembles, along with fits to the plateau regions. 



t(l.u.) 


FIG. 25: Ek2P for the lowest lying continuum ISq NN state near k = Att/L in the L = 32 ensemble, along 
with the fit to the plateau region. 



t(l.u.) 


EIG. 26: Ek2P for the continuum ISq NN states with d = (0, 0,1) near k = 0 in the L = 32 ensemble, 
along with the fit to the plateau region. 


channels. A zero of the phase shift near k ^ ^ 450 MeV is evident and occurs quite close to the 

zero of the phase shift in nature. However for k < 100 MeV, the ^Sq phase shift at rrij^ ^ 450 MeV 
and in nature become significantly different. In Fig. a linear fit is shown to the three results 
below the start of the t-channel cut, with the extracted correlated constraints on the scattering 
parameters also shown. The inverse scattering length and effective range determined from the fit 
region in Fig. |^are 

= 0.021+g|[“j , 

tm"' . = 2^96^5(5) ' («) 

The allowed region of scattering parameters is shown in Fig. and is close to containing the 
experimentally determined scattering length and effective range. Since the quark masses are un¬ 
physical, the physical values need not be contained in this region and it is interesting how close 
the current results are to those in nature. As in the ^Si-^Di coupled-channels system analyzed in 
the previous section, there is a potential self-consistency issue raised by the region of the extracted 
effective range. 
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TABLE X: Scattering information in the channel. The uncertainties are highly correlated, as can be 
seen from Fig. 


Ensemble 

iPtotl (l.U.) 

/c* cot (^(^*50) (degrees) 

All 

0 

,-n o74+(19)(26) „ „„. + (19)(26) 

<u.zm(.20)(44) ^•^‘^-(20)(44) 

243 X 64 

0 

n Qt;4+(08)(18) n+(2-0)(10.0) s+(3-0)(6-5) 

'^•3’'J^-(08)(19) '^•'^-(1.1)(1.8) ■‘■'^•®-(3.0)(6.7) 

32^ X 96 

32^ X 96 

323 gg 

0 

0 

1 

n fiQl+(09)(16) , 7+(0.5)(1.1) „„ „+(4.2)(7.0) 

■^•‘-(0.3)(0.5) ^^•'^-(4.2)(7.2) 

l-079+[S(S -3.3 +!o:SS;3 -18.3(2.6)(5.2) 
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FIG. 27: Scattering in the channel. The left panel shows fc* cot is a function of 

while the right panel shows the phase shift as a function of momentum in MeV. The thick (thin) region of 
each result correspond to the statistical uncertainty (statistical and systematic uncertainties combined in 
quadrature). The black circle in the right panel corresponds to the known bound-state result from Levinson’s 
theorem, while the dashed-gray curve corresponds to the phase shift extracted from the Nijmegen partial- 
wave analysis of experimental data m- 


VI. OBSERVATIONS ABOUT NUCLEON-NUCLEON EFFECTIVE FIELD THEORY 

ANALYSES 

A modern approach to low-energy nuclear physics rests upon the chiral nuclear forces arising 
from a non-trivial extension of yPT into the multi-nucleon sector (see, for instance. Refs. [62- 
[ 6 l]b Because of the small scales in the two nucleon systems ( 7 ^ ^ 45 MeV and l^nnl ^ 8 MeV), 
the NNEFTs are more complicated than a simple expansion in quark masses and momenta that 
defines xPT, and there are additional dynamics that must be considered. Following the initial 
developments by Weinberg [651 - I67] . much effort has gone in to understanding the construction and 
behavior of these theories. 

NNFFTs provide a powerful means with which to analyze the momentum and quark-mass 
dependences of the phase shifts and it is illuminating to consider the LQCD results presented in 
this work in their context. As is appropriate, we use KSW power counting [68H7D] in the 
channel and BBSvK power counting in], a variant of Weinberg’s power counting [65l [66] . in the 
^Si-^Di coupled channels. There are a number of reasons to undertake this investigation. The 
chiral decomposition of nuclear forces automatically requires the introduction of terms that are 
only distinguishable through variation of the quark masses. Comparison of LQCD calculations 
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FIG. 28: Scattering in the ^Sq channel below the start of the t-channel cut, < m^/4. The left panel 
shows the linear fit with the darker and lighter shaded regions associated with the statistical uncertainty and 
the statistical and systematic uncertainties combined in quadrature. The right panel shows the scattering 
parameters, and determined from fits to scattering results below the t-channel cut. The solid 

circle corresponds to the experimental values. 


at unphysical masses allows this previously unavailable “dial” to be turned in the dual expansion 
that defines chiral NNEFTs. Secondly, the full decomposition of the chiral NN forces, and thereby 
precise predictions for nuclear observables, requires knowledge of the mass dependence discussed 
above and it is essential that such calculations be performed to maximize the predictive power of 
NNEFTs. Thirdly, the current calculations enable an exploration of the convergence of NNEFTs 
with pions included as explicit degrees of freedom at relatively large pion masses. 

The quality and kinematic coverage of scattering results that have been presented is not yet 
sufficient to perform a comprehensive analysis of NNEFT matching to LQCD. Instead, we present a 
simplified discussion of the two channels to highlight some of the important features and questions 
that will need to be addressed in order to accomplish a reliable determination of the chiral nuclear 
forces from LQCD. Related discussions in the context of pionless EFTs for multi-nucleon systems 
can be found in Ref. [721I73] and implicitly in the presentation of the effective range expansion 
above. 


A. KSW Analysis of the ^Sq Channel 

The KSW power counting [68] - [70] provides a rigorous framework with which to perturbatively 
expand the two-nucleon scattering amplitude in the channel in the two small-expansion param¬ 
eters, nominally p/Aww and tyi^/Ann. Here Aww = is the natural scale of validity 

of the NNEFT. At the physical point Aww ^ 289 MeV, while at a pion mass of 450 MeV it is 
Ann ^ 350 MeV. These scales should be compared with the start of the t-channel cut from the 
next lightest meson, mp/2 ^ 385 MeV at the physical point, and mp/2 ^ 443 MeV at a pion mass 
of 450 MeV. This power counting treats the zero-derivative two-nucleon operator nonpertubatively, 
and was developed in order to correctly dehne a theory that is hnite and renormalization group 
invariant at each order in the expansion. An analysis of NN interactions at the physical point 
has been carried out to NNLO in the KSW expansion |M1 ESI ClHZSl, and we have performed 
the analogous analysis of the present LQCD results. The LO, NLO and NNLO amplitudes in 
the ^Sq channel can be found in Refs. UHl ESI dl ES], along with the relevant expansion of the 
phase shift. At LO, there is only one ht parameter, constrained by the location of the dineutron 
pole. At NLO, there are nominally two additional ht parameters, but requiring the dineutron pole 
remains unchanged reduces the number to one, ^i, while the other, ^2, can be directly related to 
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FIG. 29: The left panel shows the LQCD scattering phase shift along with the KSW NNEFT fits 
at LO, NLO and NNLO. At LO there is one parameter that is fit to recover the dineutron pole, giving 
the red-shaded region, at NLO there is one additional fit parameter, giving the blue-shaded region, and at 
NNLO there is a further fit parameter, giving the green-shaded region. The darker (lighter) shaded regions 
correspond to the statistical (statistical and systematic uncertainties combined in quadrature). The right 
panel is a scatter plot of the central values of the extracted NNLO fit parameters, ^ 1^4 over the 1 -cr range of 
the dineutron pole. The red (orange) shaded regions correspond to the statistical (statistical and systematic 
uncertainties combined in quadrature). 


Finally, at NNLO there are three more parameters, but only one parameter, ^ 4 , is independent for 
similar reasons as at NLO. Therefore, there are only three fit parameters for a complete analysis 
at NNLO. Results of fitting the LO, NLO and NNLO phase shifts are shown in Fig. The phase 
shifts at all momenta are utilized in the fits (a more complete analysis would consider the effects 
of truncations). 

Fitting the location of the dineutron bound state, the LO fit is clearly inconsistent with the 
phase shifts at higher energies, as is also seen in fits at the physical point. At NLO the fit is quite 
reasonable at the energies near the zero of the phase shift, but becomes somewhat deficient at lower 
energies. The NNLO fit is found to move closer to the LQCD results. It appears that the KSW 
expansion is converging to the LQCD results, but fits beyond NNLO are required to reproduce the 
LQCD results with an acceptable goodness-of-fit. The values of (^ 1 ^ 4 , are both of natural size, as 
can be seen in Fig. 

The resulting scattering parameters at NLO and NNLO are 

2.62(07)(16) fm ^nlo ~ 1.320(18)(38)) fm 
2.99(07)(15) fm r^o = 1-611(42)(83)) fm , (20) 


“nlo ~ 

^NNLO ~ 


From the differences between orders, it is clear that the systematic uncertainty introduced by the 
KSW expansion exceeds the uncertainties of the LQCD calculations, and orders beyond NNLO 
are required to render the “theory error” (from truncating the KSW expansion) small compared 
with the uncertainties of the calculation. As the KSW expansion is a double expansion in both 
momentum and the pion mass, the threshold scattering parameters have chiral expansions order- 
by-order in the expansion. The values of the scattering parameters extracted from fitting the KSW 
expressions differ from those obtained by fitting a truncated ERE to the phase shifts at the lowest 
two momenta and the dineutron pole, i.e. they do not lie in the region presented in Fig. This 


may indicate that the KSW expansion should not be applied to the phase shifts over the full range 
of momenta; indeed the largest two momenta have Ann- However, removing these points does 
not change the fit qualitatively due to the relative size of the uncertainties. These results could 
also indicate that the pion mass is simply too large, as it exceeds Ann- However, it does appear 
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that the expansion is converging, albeit slowly, to the calculated phase shifts. 

As the calculations have been performed at only one pion mass (previous phase shift calculations 
at = 806 MeV [HI [IT] are expected to be beyond the range of applicability of the NNEFT), 
it is not possible to isolate the explicit short-distance pion-mass dependence in (^ 1 ^ 4 , which both 
receive contributions from pion-mass independent and pion-mass dependent terms. Hence, a chiral 
extrapolation to the physical point is not feasible from this work alone. Calculations that are 
currently underway will provide results at a lower pion mass, from which predictions at the physical 
point will become possible. 


B. BBSvK Analysis of the Coupled Channel 

BBSvK power counting m is similar to Weinberg’s power counting [Ml [Ml, and is an appro¬ 
priate scheme to use in the case of the ^Si-^Di coupled channels. An NN interaction (two-particle 
irreducible) is derived using the familiar rules of yPT and, due to the infrared behavior of the two- 
nucleon system, is iterated to all orders with the Schrodinger equation to generate the bound-state 
pole(s) and scattering amplitude(s). See Ref. [77] for a review. 

At LO in Weinberg’s power counting, the NN interactions are determined by momentum- 
independent and quark-mass-independent two-nucleon contact interactions and by one-pion ex¬ 
change (OPE). However, the short-distance nature of the tensor force, resulting from OPE, gener¬ 
ates renormalization-scale dependence in the D-waves that requires the presence of a counter term 
at LO, and BBSvK is the simplest power counting to remedy this situation. At NLO in the counting, 
there are contributions from pion-loop diagrams, from momentum-dependent two-nucleon contact 
interactions, and from insertions of the light-quark mass matrix into momentum-independent two- 
nucleon contact interactions. With the parameters in the meson sector, e.g. /tt, fixed to the 
results of other LQCD calculations at similar quark masses, there is one free parameter at LO in 
BBSvK counting - the coefficient of the momentum-independent two-nucleon contact interaction. 
This is common to both the S-waves and D-waves. 



0 100 200 300 400 500 

k (MeV) 


FIG. 30: The coupled channels scattering phase shift, 6 ia^ along with BBSvK fits at LO and NLO*. 

At LO there is one parameter that is fit to recover the deuteron pole, giving the red-shaded region, while 
at NLO* there is one additional fit parameter, giving the blue-shaded region. The darker (lighter) shaded 
regions correspond to the statistical (statistical and systematic uncertainties combined in quadrature). 


At NLO, the expansion becomes more complicated with different interactions in the S-waves and 
D-waves. Without being able to separately resolve the 6 ia and 61 ^ phase shifts, only the common 
terms can be determined. To this end, we have defined NLO* to be LO with the inclusion of the 
leading momentum-dependent two-nucleon contact interaction that is also common to both the 
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S-waves and D-waves, but omitting other NLO contributions. NLO* introduces a single additional 
parameter beyond LO. 

The results of fitting the LO and NLO* parameters to the results of our LQCD calculations 
are shown in Fig. 30 ^ The LO fit to the deuteron binding energy leads to phase shifts that 


significantly over estimate the LQCD results (this is slso seen in analyses at the physical point). 
However, by including the contact-p^ interaction, relatively good agreement is found in the NLO* 
fit to all the LQCD phase-shift extractions, with the exception of the lowest energy point. The 
values of the scattering parameters resulting from the fits are 


^LO ~ 
^NLO* ~ 


1.94(09)(17) fm 
2.72(22)(27) fm 


= 0.674(17)(29)) fm 
~ 1.43(12)(13)) fm 


( 21 ) 


which are consistent, within uncertainties, with those obtained in the ^Sq channel with KSW 
counting. It is interesting to note that the ratio of scattering length to effective range is a/r ^ 2, 
as was found to be the case at the SU(3) symmetric point [HllTT] . 

A feature of BBSvK counting is that predictions can be made for the mixing parameter, ei 

and or Ii and in the more familiar Stapp [78] parameterization of the S-matrix. These 

are shown in Fig. and it is important to keep in mind that the coefficients determined from 
the deuteron pole and S-wave phase shift, contribute to both these quantities. While the D-wave 



FIG. 31: The left panel shows the ^Di scattering phase shift, , in the Stapp Parameterization m 
along with BBSvK fits at LO and NLO*, while the right panel shows the mixing parameter, ei. The darker 
(lighter) shaded regions correspond to the statistical (statistical and systematic uncertainties combined in 
quadrature). 


phase shift is only slightly modified by the NLO* interaction, 6i is changed dramatically. In this 
initial investigation, the range of the square well interaction has not been varied and estimates of 
contributions from higher orders have not been included. It is clear that the “theory error” due to 
truncation of the BBSvK expansion is large for ei, but not for the D-wave phase shift. In fact, this 
expansion of 6i is found to be less convergent at this pion mass than at the physical point m- 


^ A square-well with a radius of = 0.30 fm has been used to regulate the interaction at short distances. Previous 
work HU shows that the observables have corrections that depend only on positive powers of R (after refitting 
coefficients), as expected from a Wilsonian renormalization group analysis in the limit R ^ 0. 
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VII. CONCLUSIONS 

Recovering the experimentally known properties of the two-nucleon systems, such as the 
deuteron bound state, the dineutron virtual-bound state and scattering observables, from QCD 
represents a major challenge for Lattice QCD calculations. Once verified by comparison to known 
experimental extractions, LQCD calculations hold the promise of refining our knowledge of these 
systems beyond what is possible experimentally, particularly in the neutron-neutron system and 
more exotic processes involving hyperons. LQCD calculations have steadily developed in recent 
years and in the near future calculations of multi-nucleon systems with physical quark masses will 
be available m- Eventually these calculations will also include the effects of isospin-breaking 
and QED. In this work, we report the results of calculations of nucleon-nucleon interactions in 
the ^Si-^Di coupled channels and the ^Sq channel at a pion mass of rrij^ ^ 450 MeV in three lat¬ 
tice volumes and at a single lattice spacing. ^ Both the deuteron and dineutron are found to be 
bound at this pion mass, consistent with expectations based upon previous calculations. The phase 
shifts in both channels are determined at a few discrete momenta and, in both channels, a zero 
in the phase shift is found to occur near the momentum at which a zero is observed in nature. 
Calculations of increased precision and kinematic coverage will further our understanding of the 
two-nucleon systems at this set of quark masses. Further calculations at other quark masses will 
enable direct comparison with experimental extractions and will elucidate important features of 
the chiral nuclear forces that are not accessible in experiment alone. 
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